One-equation k RANS: wall condition, closure fixes, terrain wall distance, Canonical_RANS regression cases - #3971

…e step

With erf.dirichlet_k the surface layer writes the AL01 Eq. 16 value of k
into the first cell of S_old once per step. The logical BC for RhoKE was
then switched to ext_dir, so the diffusion applied a one-sided stencil
against a zero ghost cell and drained a fifth of the wall k before the
step ended; the converged wall k sat 20 % below the AL01 value and the
implicit KE solve was silently switched off by the BC sanity check.

Keep foextrap for RhoKE at a surface-layer wall (ghost cell equals the
first cell, zero flux through the wall face), restore the first-cell
value from S_old after every RK stage in erf_slow_rhs_post, and pin the
bottom row of ImplicitDiffForStateLU_{N,S,T} for RhoKE so the implicit
solve treats it as a Dirichlet neighbour.

Neutral_ABL_Flat, 12 h: KE(0)/u*^2 = 3.232 against 3.2325 (was 2.60),
k_start/k_end = 1.000 (was 1.244), interior profile unchanged, restart
bit-exact, and the wall-cell check now runs in the CTest smoke mode. The
implicit path was exercised with a compressible variant since anelastic
disables it. erf.dirichlet_k keeps its default of false.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

Without the Dirichlet wall value the first cell cannot resolve the
near-wall shear production and k there settles at about half the
Axell & Liungman equilibrium u*^2/Cmu0^2 (12 h neutral run: 1.59 against
3.23), while the mean wind still follows the log law because MOST
supplies the stress. Warn at startup in that configuration; the default
stays false. Every deck that enables kEqn now sets the flag (Askervein
added; the canonical deck already had it).

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

…table length, unit tests

Factor the Axell & Liungman closure relations into ERF_RANSClosure.H
(namespace AL01) so ComputeTurbulentViscosityRANS and the unit tests
evaluate the same expressions; the refactor is bit-identical to phase 2.

Validate Cmu0, Cb, max_geom_lscale, Rt_crit <= 0, Rt_min < Rt_crit and
Rt_min > -3.6 (the poles of Eqs. 31-32) at input time. Add
erf.tke_floor (runtime floor on k; distinct from erf.tke_min, the initial
value) and the opt-in erf.implicit_tke_dissipation, which drops the
explicit sink from the source and divides the update by 1 + dt c with
c = diss_old / (rho k)_old, half-weighted on the anelastic stage 1.

Two robustness changes came out of the new gtests and the review against
the paper and the Kynema KLAxell implementation: the Burchard & Petersen
smoothing is rewritten as Rt_crit + a x / (x + a), the same function
without the cancellation that returned -2 for Rt = -1e16 and +1.4e14 for
Rt = -1e30; and the unstable length evaluates Eq. 28 once with the
smoothed Rt from the geometric length, bounded by about 1.31 l_g, instead
of the two-pass corrector that iterated the fixed-point map of Eq. 26,
which has no fixed point in strong convection (AL01 p. 78).

Also: fill the derived diffusivity components on the tilebox only (the
wrapper fills ghost cells), drop unused parameters, correct the surface
heat-flux index in the TKE-source comment, and use the level's own
periodicity in the wall-distance Poisson solve.

Neutral_ABL_Flat: all checks unchanged to the digits shown; the 4x-dt
run showed the anelastic integrator's explicit vertical diffusion is the
real time-step limit (dz^2/2K about 19 s), not the dissipation; a
compressible dt = 60 s pair with implicit vertical diffusion passes with
the dissipation explicit and implicit.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

…er anelastic

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

…ap, stable and convective decks

Two opt-in closure options: erf.rans_consistent_diffusivities makes the
horizontal heat and every scalar and moisture diffusivity follow the AL01
scalar stability function (rho cmu' sqrt(k) L) like the vertical heat
diffusivity, instead of Pr_t and Sc_t times the eddy viscosity; and
erf.rans_lscale_from_pblh caps the geometric length at kappa 0.1 zi from
the surface layer's MYNN25 PBL-height diagnostic, clamped to
[erf.rans_lscale_min, erf.max_geom_lscale]. The neutral deck is
bit-identical with both off.

Two new Canonical_RANS cases with numeric check scripts and CTest smoke
entries: Stable_ABL_Flat (GABLS1, 9 h) passes every check first time,
u* 0.244, jet 1.23 Ug at 154 m, depth 134 m; Convective_ABL_Flat
(MS94-B sounding, 0.24 K m/s, 4 h) closes the column heat budget to
0.02 %, keeps the inversion at 1020 m and carries the AL01 buoyancy term
at the wall. It runs at dt = 2 s: at 5 s the eddy viscosity reaches the
explicit vertical-diffusion limit of the anelastic integrator (40 m2/s)
and the run aborts, the same limit phases 9 and 10 address. The
mixed-layer spread tolerance (2 K) reflects the superadiabatic lapse a
local-K closure keeps without countergradient transport; Kynema's
KLAxell has none either and its Prandtl function matches cmu'/cmu.

Shared check code moves to rans_checks.py; the dissipation-lag check
allows 10 % in smoke mode (early transient) and 5 % in physics mode.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

… deck

The Poisson wall distance (Tucker 2003) took the cell's low-face fluxes
as its gradient, half a cell off centre in every direction, which
overstates |grad phi| by dz/2 and shortens every distance by z dz/(2H) on
a flat mesh, 0.78 % for 64 cells. Form a cell-centred gradient instead
from centred differences with the cell-centre terrain metrics (chain rule
for a mesh deformed in z); it needs one ghost cell of phi and the cell's
own nodes. On a flat fitted mesh the distance is now exact to 1e-6 m
above the first cell (1.5 cm there, the odd-reflection Dirichlet ghost);
on the ridge the mean error drops from 1.5 % to 1.0 %.

New Canonical_RANS case Neutral_Hill_2D: periodic Witch of Agnesi ridge
(h 100 m, L 500 m) on a terrain-fitted mesh under MOST, one cell wide in
y, anelastic with the MLMG projection. Its check script compares the wall
distance with the exact distance to the curve, and after 6 h measures the
crest speed-up (0.55, 0.45, 0.39 in the lowest cells against the
Jackson-Hunt estimate 0.4) and the upstream log law (within 4 %). A
flat-fitted variant of the same deck (prob.hmax = 1e-6) is the second
CTest entry; both pass. erf_plotfile.py gains a full-field reader for the
terrain checks.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

…estart on terrain

Add erf.wall_dist_type = terrain_height: the height above the local
surface projected on its normal, no linear solve (after the amr-wind
immersed terrain and Kynema). It is exact to 1e-10 on a flat fitted mesh
and closer to the true distance than the Poisson solve on both hills
(mean 0.02 % and 0.01 % against 1.0 % and 0.3 %); the hill decks use it
and three _Poisson CTest variants keep the Poisson path covered. The
Poisson solve is now posed in the positive-definite form MLABecLaplacian
documents (same iterates as before).

New case Neutral_Hill_3D: a periodic radial Witch of Agnesi hill at
dx = dy = dz = 40 m; after 4 h the crest speed-up is 0.41, 0.28, 0.22 in
the lowest cells against the axisymmetric estimate 1.6 h/L = 0.32, and
the upstream log law holds within 10 %. Restart is bit-exact on the 2D
and 3D terrain decks; the Askervein deck runs clean.

Recorded in the plan, outside its scope: on a 3D fitted mesh with dz
different from dx the initial field's divergence is a deterministic
1.788e139 before the first projection and the run aborts, independent of
the closure, the boundaries, the terrain source, stretching, the box
layout and the projection solver, while Askervein runs at the same aspect
ratios; a separate uninitialised read in the w boundary fill trips the
invalid-operation trap in the first advance. Reproducer in PLAN.md.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

…agnostics

Add Docs/sphinx_doc/theory/RANS.rst, describing the Axell & Liungman
one-equation closure as implemented: the eddy viscosity and diffusivity,
the k equation and its dissipation, the geometric and stratified length
scales with the bounded unstable form, the stability functions with the
cancellation-free smoothing, the wall condition on k, the two wall
distances on a fitted mesh, the limitations and an input table. Link it
in the THEORY toctree, document the Canonical RANS regression cases in
RegressionTests.rst, and add a top-level README for the case directory.

Three new EddyDiff components, written by the closure, expose the
smoothed turbulent Richardson number and the two stability functions as
the plot variables Rt, cmu and cmu_prime; the neutral smoke check
verifies them against AL01 Eqs. 31-32 to round-off and Kmv against
rho cmu sqrt(k) L.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

…istency check

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

# Conflicts:
#	Tests/CTestList.cmake

…ain cases on FFT

No CI configuration builds with ERF_ENABLE_FFT, and every Canonical_RANS
deck ran the anelastic projection with the FFT solver, so all nine RANS
entries aborted at start-up in CI (the 2026-09-10 run: 9 of 9 on macOS,
Linux GCC and Windows). The flat decks abort in ReadParameters on
erf.use_fft; the 2D hill deck sets use_fft = false but its general-terrain
projection is the FFT-preconditioned GMRES, which has no non-FFT path.

- The three flat entries pass erf.use_fft=false (MLMG projection) and run
  in every build; the six terrain-fitted entries are registered only under
  if(ERF_ENABLE_FFT).
- The test command echoes the tail of the run log when the executable
  exits non-zero, so ctest -VV shows the reason instead of nothing.
- RegressionTests.rst and the Canonical_RANS README say so.

Checked with a Debug build without FFT (3 entries registered, 3 pass), a
configure with MPI off (3 entries), and a Release FFT build (9 pass).

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

…; correct the dz != dx notes

The terrain_height wall distance read z_nd(:,:,klo) from every box, but
boxes above a z split of the BoxArray do not hold the surface nodes, so
their distance was garbage (778 m error on the 3D hill deck at
amr.max_grid_size_z = 10, clamped to epsilon). Gather the surface slab onto
every box with a ParallelCopy first, and guard the gathered slab with a
reduction over its own ghost width (a scalar-ghost MultiFab::max would read
past the slab, which has no z ghosts). With the fix the split and unsplit
layouts give bit-identical wall distances, and the unsplit distance is
unchanged.

The 1.788e139 pre-projection divergence recorded in PLAN.md phase 6, the
3D hill README, RESULTS.md, the deck and the Poisson solver comment as a
dz != dx failure was the same z-split layout (unfilled momenta ghost faces
in the initial projection and duplicated planar surface-layer arrays,
fixed in erf-model#3970), not the aspect ratio; the text now says so.
The Poisson wall-distance multigrid does still diverge at dx = 2 dz.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>